Model Flames in the Boussinesq Limit: The Effects of Feedback 
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We have studied the fully nonlinear behavior of pre-mixed flames in a gravitationally stratified 
medium, subject to the Boussinesq approximation. Key results include the establishment of criterion 
for when such flames propagate as simple planar flames; elucidation of scaling laws for the effective 
flame speed; and a study of the stability properties of these flames. The simplicity of some of our 
scaling results suggests that analytical work may further advance our understandings of buoyant 
flames. 
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I. INTRODUCTION 

In several areas of research, the feedback of a propa- 
gating diffusive (pre-mixed combustion) flame on a fluid, 
and the consequent effects of the flame itself, is of consid- 
erable interest. In the astrophysical context, for exam- 
ple, the speedup of nuclear reaction fronts of this type 
in the interior of white dwarf stars is thought to be one 
possible way that such stars undergo thermonuclear dis- 
ruption, e.g., a Type la supernova (cf. 0, 0, 0, El IE El ) • 
Much of the literature on this subject has focused on the 
speedup of such flames for prescribed flows, and substan- 
tial advances have been made in this regard recently ; 
this is the "kinematic" problem, in which one seeks to 
establish rigorous limits on flame speedup in the case in 
which there is no feedback onto the flow. The aim of this 
paper is to study the simplest case of feedback, namely 
that which occurs when a flame propagates vertically, 
against the direction of gravity. As described extensively 
in the previously cited literature, it is generally believed 
that under such circumstances, the flame front is likely 
to become distorted by the action of the Rayleigh- Taylor 
instability, and thus achieves speedup; these calculations 
have been largely illustrative, and based upon simula- 
tions using fully-compressible fluid dynamics (e.g., @) 
and fairly realistic nuclear reaction networks. 

Here, we focus on a much simpler problem; we study 
such flames in the Boussinesq limit (leading to a far sim- 
pler computational problem) and for highly simplified 
reaction terms (avoiding the complexities of realistic nu- 
clear reaction networks). In this way, we are able to 
isolate the various effects which lead to flame speedup, 
which is particularly important if one is to connect such 
simulations to the extant analytical work on this subject 
( e -g-> 0E3 E3) Indeed, an important motivation for this 
work is to elucidate simple scaling laws — if they exist 
— in order to suggest further analytical studies. 

Our paper is structured as follows: In the next sec- 



tion, we describe the specific physical problem we wish to 
study, establish the equations to be solved, and describe 
the method of solution. In §111, we present our results, 
and in §IV, we provide a summary and discussion. 



II. THE PROBLEM 

The effect of gravity on the temperature distribution 
in a reacting incompressible fluid with thermal diffusivity 
k, viscosity /i, and density p can be described by the set 
of Navier-Stokes and advection-diffusion-reaction equa- 
tions, 



<9v 

dT 
~dt 



Vp + /iV 2 v + pg, 
R(T), 



kV 2 T 
0. 



(1) 



where v is the fluid velocity and, without loss of gener- 
ality, the temperature T has been normalized to satisfy 
< T < 1. The thermal diffusivity and viscosity are as- 
sumed to be temperature- independent, and density vari- 
ations are assumed to be small enough to be described by 
the Boussinesq model, e.g., p(T) = p a + (Ap/p )T. The 
model Q can be derived from a more complete system 
under the assumption of unity Lewis number Le (the ra- 
tio of thermal and material diffusivities), and this article 
addresses only the Le = 1 case. 

The Boussinesq model is the simplest system exibit- 
ing buoyancy effects (and thus allowing for feedback to 
the flame) without introducing the complexities associ- 
ated with the presence of sound wave and stratification 
of background "atmosphere" . Because our intentions is 
to elucidate basic priciples, rather than realistically mod- 
eling specific physical situations, we view our approach 
as sufficient for the chosen task. 

We consider a reaction term of Kolmogorov-Petrovskii- 
Piskunov (KPP) type ^lj, °f the form 



R(T) = aT(l-T)/4. 



(2) 
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where a is the (laminar) reaction rate. This reaction 
form has an unstable fixed point at T = 0, the "un- 
burned" state, and a stable one at T = 1, the "burned" 
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initial velocity v =0 



FIG. 1: A typical initial state of a flame calculation. 



state. Thus a fluid element with positive temperature will 
inevitably evolve to the burned state in a characteristic 
time of order l/a. As is well-known from the combus- 
tion literature, the temperature equation from the sys- 
tem above admits — for a stationary fluid, and in the 
absence of gravity — one-dimensional solutions in the 
form of burning fronts propagating with laminar burning 
speed s , and with characteristic flame thickness 5, 



(3) 



If it is further assumed that T — ► 1 as y — > — oo, and 
T — ► as y — > +oo, then the front propagation is in the 
positive y direction. 

It is convenient to adopt the front thickness S and the 
inverse reaction rate a -1 as the units of distance and time 
respectively. In these units the problem control parame- 
ters are the Prandtl number Pr and the non-dimensional 
gravity G, 



Pr=^ 

K 



(4) 



where v is a kinematic viscosity v — fi/ p . In addi- 
tion, the system is characterized by a number of length 
scales specifying the initial state, which are in our case 
the dimcnsionless amplitude A and the dimensionlcss 
wavelength L of the initial flame front perturbation, 
f(x) = acos{2-Kx/l), 



A = a/5, 



L = l/6. 



(5) 



The vertical size of the computational domain was kept 
large so as to avoid effects due to the upper and lower 
walls of the computational box; in all cases, we have ver- 
ified that such artifacts are not present. For this reason, 
the box height does not enter as a problem parameter. 
The initial velocities are set to zero, and most computa- 
tions were carried out for Pr = 1. A typical initial state 
of our flame calculation is shown in Fig. ^ 



Because we focus on the two-dimensional problem, it 
is convenient to re- write Eqns. in the stream function 
and vorticity formulation in dimcnsionless form, 



dui 
~dt 



= —v • Vw 



PrV^w 



G 



dT 

dx 



dT 
~dt 



= -v • VT 



V 2 T 



-T(l-T). 



(6a) 



(6b) 



using S and 8 / s Q as units of length and time respectively. 
Here v is the non-dimensional velocity and ui is the non- 
dimensional vorticity (w e V x v = V 2 i/i). We solve the 
system (Eqns. numerically. The solution is advanced 
in time as follows: a third order Adams-Bashforth inte- 
gration in time advances to and T, where spatial deriva- 
tives of lu and T are approximated by fourth-order (ex- 
plicit) finite differences. The subsequent elliptic equation 
for ip is then solved by the bi-conjugate gradient method 
with stabilization, using the AZTEC library [l2|]. Finally 
we take derivatives of ip to update v. 

The resolution of the simulations was chosen to fully re- 
solve the laminar flame structure. For the KPP reaction 
term 0, the laminar flame thickness is approximately 
12 5, and the grid spacing Ax = Ay = 1 (in the units of 
5) was used in most of the computations. The laminar 
flame speed computed at this resolution agrees with the 
theoretical value to within 1%. This corresponds to at 
least 16 zones per wavelength (of the initial perturbation) 
which is sufficient to resolve the flow. Most of the compu- 
tations were executed on a domain half the width of the 
initial perturbation wavelength, and reflecting boundary 
conditions were applied. 

Simulation times of t — 200 - 500 (in units of S/v ) 
were required to measure the bulk burning rate, on com- 
putational grids ranging from 8 x 256 for L = 16 to 
64 x 2048 for L — 128. Larger domains were necessary 
to obtain velocity fields (e.g. 64 x 3072 for L = 128), 
to avoid the influence of boundary conditions at the top 
and bottom. Fortunately, only the velocity in the re- 
action region affects the shape of the flame front and, 
consequently, the bulk burning rate, so slight errors in 
estimating velocities far away from the front due to up- 
per and lower boundaries do not affect our results. The 
comparison with linear analysis was done using the same 
resolution and domain sizes up to 512 x 4096. 

By its nature, this study was comprised of a large num- 
ber (~ 250) of simulations each representing a data point, 
as opposed to a close examination of just a few simula- 
tions as in a case-study approach. Confidence in the nu- 
merical accuracy was gained at the cost of a small num- 
ber of additional test simulations. Some simulations were 
repeated using lower and higher resolutions, domains of 
different sizes in vertical direction, and with several wave- 
lengths across the width of the domain. Special attention 
was devoted to simulations with different Prandtl num- 
bers, to ensure that both diffusive and viscous scales were 
resolved. 
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Finally, a comment regarding the two-dimensional na- 
ture of our simulations. In a recent study of the closely re- 
lated Rayleigh- Taylor instability ^| , Young et al. specif- 
ically compared the behavior of fingering and mixing in 
two- and three-dimensional flows, with the result that 
while the specifics, e.g., finger growth rates, were quite 
sensitive to dimensionality, the phenomenology neverthe- 
less turned out to be rather similar. Flames do however 
introduce a very useful physical simplification into the 
Rayleigh- Taylor problem: Because flames consume all 
density features at the flame front with scales smaller 
than the flame thickness, the Rayleigh- Taylor problem is 
"regularized" by the burning process even in the limit 
of vanishing viscosity. For this reason, a key difference 
between 2-D and 3-D - namely, the difference between 
small-scale turbulent structures in two and three dimen- 
sions - is sharply reduced in the burning case. The re- 
maining difference between 2-D and 3-D is then mostly 
related to the difference in propagation speed between 
buoyant parallel rolls (the 2-D case) and buoyant tori 
(the 3-D case), with tori propagation more quickly, i.e., 
we would expect 3-D flames to propagate more quickly 
than 2-D flames, all other things being equal. We plan 
to explore this point in future three-dimensional studies 
of flame propagation. 



III. RESULTS 

In this section, we discuss the results of our calcula- 
tions, focusing successively on the bulk burning rate, the 
evolution of the burning travelling front, and the ultimate 
transition to a travelling (burning) wave. Our central in- 
terest is in disentangling the dependence of the flame 
behavior on the key control parameters of the problem. 



A. Travelling wave flame 

For a wide range of parameters, we were able to con- 
struct a sufficiently large computational domain that we 
could observe travelling waves of the temperature dis- 
tribution, propagating with constant speed. Depending 
on simulation parameters, the initial perturbation either 
damps (e.g., the flame front flattens) or forms a curved 
front. The flat front moves in the motion-free (in the 
Boussinesq limit) fluid, has laminar front structure, and 
propagates with the laminar front speed. 

The typical curved front is shown in Fig. it has the 
wavelength of the initial perturbation and is character- 
ized by narrow dips (lower apexes), where the cold fluid 
falls into the hot region, and by wide tips (upper apexes), 
where buoyant hot fluid rises into the cold fluid. In the 
initial stages, the evolution pattern is similar to bubble 
and spike formation during the Rayleigh- Taylor instabil- 
ity H3> in later stages, small scale structures are 
consumed by the flame, and, finally, the flame evolves 
toward the travelling wave solution shown in Fig. [21 The 



L = 32 L=128 
G = 1 G = 1/4 




FIG. 2: Travelling wave isotherms (T = 0.1 and T = 0.9) 
and streamlines for two system with different simulation pa- 
rameters. Note that the system on the right has been rescaled 
by a factor of 1/4 both horizontally and vertically. 



shape of the stable front is determined by gravity, G, and 
wavelength, L, and can be characterized by two verti- 
cal length scales associated with the spatial temperature 
variation (hr) and the spatial velocity variation (hy) of 
the flame. The speed of the curved front is always higher 
than the laminar flame speed, because of the increase in 
the flame front area and transport. Finally, we notice 
that the streamlines in Fig. |2 indicate that the flow un- 
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FIG. 3: Bulk burning rate (travelling wave speed) s as func- 
tion of wavelength L for different values of gravity G. 

derlying the propagating flame is characterized by rolls 
propagating upward. 

One of our primary interests is to quantify the effects of 
variations in wavelength and gravity on the flame speed. 
It is convenient to define the speed of the travelling wave 
flame by the bulk burning rate 0, 

1 f l dT(x,y,t) 
S{t) = lh^^ dxdy] (7) 

this definition has the considerable advantage that it re- 
duces to the standard definition of the flame speed when 
the flame is well-defined, and it is accurate to measure 
even for cases where the burning front itself is not well- 
defined. Henceforth we refer to it simply as the flame 
speed. 

Our first result (shown in Fig. [3J| is that the flame 
speed increases with wavelength L and with the gravi- 
tational acceleration G, and is independent of the initial 
perturbation amplitude A. More specifically, the flame 
becomes planar and moves at the laminar speed (s = s a ) 
if G is smaller than some critical value G CI ; if G lies 
above this critical value, the flame speed can be fit by 
the expression, 

s = s y/1 + ki(G - d)L , (8) 

where k\ » 0.0486 is obtained from measurements de- 
rived from the simulation data. The second tuning pa- 
rameter, Gi, was found to be a function of the pertur- 
bation wavelength (Fig. gj), Gi = 8(2n/L) 1 - 72 . For a 
relatively wide range of parameters, Eq. JSJ describes ex- 
perimental data well, but must be applied with caution 
near the cusp at G = G\ shown in Fig. |3J Roughly 
speaking, this cusp can be interpreted as the transition 
between the planar and curved flame regimes, Gi ~ G CT ] 
closer investigation of the transition region shows that 
G cr < Gi , and that the fit (Eq. |8J) underestimates the 
flame speed in this transition region (Fig. EJ. 




2n/L 

FIG. 4: Transitional point Gi as a function of wavelength. 

The behavior near the transition is discussed in detail 
in the theoretical work carried out by Berestycki, Kamin 
& Sivashinsky |l0j | . They derive the one-dimensional evo- 
lution equation for the front interface y(x) and prove 
mathematically the following properties of y(x), relevant 
to our case: (1) the existence of G cr ~ (2w/L) 2 such that 
there is no nontrivial solution for G < G cr (i.e. the front 
is flat for G < G cr ); (2) the existence of G cr * = 4G cr such 
that for G > G cr * there are two symmetrical (curved) so- 
lutions y + (x) and y~(x) which are stable, and any other 
solution including the trivial is unstable; (3) metastabil- 
ity of any solution except y + (x) and y~{x) in the range 
G cr < G < G cr *, and convergence of this metastable 
solution to either y + (x) or y~{x). Also, based on the 
derivation in [To| . it can be shown that the flame speed 
in the metastable regime scales as follows [T^. 

(s/s - 1) cx (G - G cr ) 2 as (G-G cr )^0. (9) 

3 | , , , 1 
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FIG. 5: Amplitude of the stable front as function of gravity 
for the wavelength L = 32. The scaling relations shown here 
are discussed in the text. 
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Our simulations confirm the dramatic increase of sta- 
bilization times close to the critical gravity value G cr . 
For this reason, it is very difficult to obtain reliable re- 
sults regarding the flame speed in this transition regime. 
Even detecting the critical point takes significant com- 
putational effort (Fig. EJ); measuring the velocity, which 
in this parameter regime differs from s by a very small 
amount, is harder still. 

However, the transition is sharper and is easier to see 
when studying the vertical distance between the upper 
and lower apexes of the flame, hx, measured by the ex- 
pression, 

/oo 
(T(0)-T(l/2))dy. (10) 
-oo 

In the limit of large wavelengths (L ^> 1), the transition 
occurs at small values of gravity, and the flame speed is 
determined by a single parameter, the product LG. If, in 
addition, the pro duc t LG is large, the flame speed scales 
as s/s 0.22^/LG. This result is in good agreement 
with the rising bubble model |17| which, in the Boussi- 
nesq limit, predicts s/s = y/LG/6-rr ~ 0.23^/LG for a 
2-D open bubble [l^. We further observe that in the 
large wavelength limit, the hx/l ratio obeys the same 
scaling (Fig.EJ. 

We note that the flame structure shares features of 
flame propagation from both shear and cellular flow. For 
instance, the temperature distribution closely resembles 
that of a flame distorted by a shear flow, while the ve- 
locity distribution resembles that inside an infinitely tall 
cell. The flame speed in the shear and cellular flow is 
determined by the flow speed and by the length scale of 
the flow (period of shear or cell size) flflj . In particu- 
lar, in both cases the flow speed scales with maximum 
flow velocity as s oc v 7 ^^, with n = 1 for burning in 
the shear flow and n = 1/4 for burning in the cellular 
flow. Similarly, we have tried to determine whether the 
flame speed relates to the maximum velocity of the flow 
when flow and flame are coupled through the Boussi- 
nesq model. The available data (shown in Fig. EJ) do not 
demonstrate a power law dependence with a single well- 
defined power. Furthermore, the dependence on L is not 
as dramatic as in the cases of shear or cellular flows. 



B. The thin front limit 

The thin front limit is particularly important for de- 
veloping models of flame behavior. For many applica- 
tions — especially in astrophysics — resolving flames 
(by direct simulation) is prohibitively expensive, and un- 
derstanding flame propagation in the limit in which the 
flame front becomes indefinitely thin (when compared to 
other length scales of the application) is critical for de- 
signing flame models. Of course, this same limit is of 
intrinsic mathematical interest. 

Particularly important is the dependence of the flame 
speed on the wavelength of the front perturbation in the 




FIG. 6: The flame speed as function of maximal flow velocity. 

thin front limit. We have already pointed out that in- 
stabilities with larger wavelengths have higher travelling 
wave speeds, so that eventually the instability with the 
largest wavelength allowed by the system dominates. (In 
our non-dimensionalization, this is the instability with 
the highest ratio of wavelength to laminar front width) . 

In this context, it is convenient to switch from 
our "laminar flame units" to the so-called "G-equation 
units". The G-equation is a model for reactive systems 
where very low thermal diffusivity is exactly balanced by 
high reaction rate (see e.g., H(|). The diffusion and reac- 
tion terms in the temperature equation are replaced by 
a term proportional to the temperature gradient, 

dT 

— +vVT = So VT , 
ot 

so that the front propagates normal to itself at the lam- 
inar flame speed s . The Boussinesq fluid model, com- 
bined with the G-equation flame model, has the following 
physical parameters: (1) flow length scale I, (2) laminar 
flame speed s , (3) gravity <?, and (4) fluid viscosity v. 
Choosing I and l/s as the length and time units, the 
governing non-dimensional parameters are g = glj s 2 and 
v = v/ls ] the corresponding parameters in the lam- 
inar flame unit system are g = LG and v = Pr/L. 
Note that in the limit L — > oo while keeping Pr = 1, 
the Navier-Stokes equation becomes the Euler equation 
and v —> 0, leaving only one parameter in the system, 
g = (Ap/ P o)gl/ S l = LG. 

In our simulations Pr = 1, so it is not surprising that 
for large L, almost all aspects of the system are well 
characterized by the LG product alone. For example, the 
formula for the bulk burning rate, s/s = VI + k\LG for 
G > G\ well describes our experimental results. Next, 
consider the travelling wave solutions shown in Fig. |3 
these two systems have the same LG product, with L = 
32 and L — 128, and move at the speed s/s a = 1.34 
and s/s — 1.51 respectively. The wavelength here is 
comparable to the laminar front thickness (indicated by 
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FIG. 7: Travelling wave solution for two systems with LG = 
32, with L = 128 (dashed lines), and with L — 256 (solid 
lines). The isotherms T = 0.1 and T — 0.9 are shown in the 
middle panel. The top panel shows the temperature profiles 
and vertical velocities (along g) at x = 0.5 (upper flame front 
apex); the bottom panel plot shows the same things at x = 
(lower flame front apex). 

the two limiting isotherms T = 0.1 and T = 0.9). Still, 
the front shape as well as flame speed and fluid velocities 
are very similar. 

One can see similarity more clearly in Fig. [7] (mid- 
dle panel), which compares systems with L = 128 and 
L = 256. The agreement between bulk burning rates 
is very good {s/s = 1.51 and s/s D = 1.57). The 
match between the two integral measures hr is weaker 
(hr/l — 0.83 and hr/l — 0.75), suggesting that the sys- 
tems in consideration are still far from the infinitely thin 
front limit, but this is apparent from the distance be- 
tween limiting isotherms. We have also compared the 
temperature and velocity profiles at the upper and lower 
apexes of the flame (Fig.0 the top and the bottom pan- 
els). The velocity is — as expected — essentially zero 
well ahead of the temperature front, but significant mo- 



TABLE I: Three simulations with LG = 32 discussed in the 
text. 



setup 


L 


G 


s/So 


hr/l 


^max / So 


(a) 


32 


1 


1.34 


0.98 


4.40 


(b) 


128 


1/4 


1.51 


0.83 


5.06 


(c) 


256 


1/8 


1.57 


0.75 


5.01 




FIG. 8: The isotherm T — 0.5 during the instability growth 
phase, shown for three systems with LG = 512 but different 
L (top: L = 64, middle: L = 128, bottom: L = 256). The 
initial amplitude is a/l = 1/8, and snapshots are taken at 
times t( So /l) = 0, 1/16, 2/16, 3/16, 4/16. 



tion extends far behind it; the absolute maximum veloc- 
ity is located in the vicinity of the lower apex and is re- 
lated to the bulk burning rate (Fig.^J. By examing the 
detailed velocity profiles we find that velocities at the 
flame front also obey the LG product scaling and, to- 
gether with the temperature distribution, determine the 
bulk burning rate. However, the velocities well behind 
the front can be quite different for two systems with the 
same LG product (cf. Fig.0). 

Finally, consider the temperature during the instabil- 
ity growth phase, shown for three different cases (with 
LG = 512) in Fig. [S] Although the wavelength to lam- 
inar front thickness ratio affects small scale features, we 
again clearly see the similarity scaling connecting these 
solutions. 

As we have shown above, the dependence on a single 
parameter, namely the LG product, in the infinitely thin 
front limit follows from dimensional analysis; and for rea- 
sonably thin fronts, we were able to confirm the LG prod- 
uct scaling. At the same time, we noticed that the length 
of the velocity variation, hy, does not scale with LG = g. 
It is reasonable to assume that hy is controlled by the 
other parameter, namely, the non-dimensional viscosity 
v = Pr/L, which is essentially zero in the thin flame 



7 



1 

0.8 

fa 0.6 

g o.4 

0.2 





-10 



y/l 




FIG. 9: Vorticity generation in the roll (solid line) and vortic- 
ity fluxes through the separatrices between the rolls (dashed 
line), as a function of height y for a flame with L — 32 and 
G — 1 (top panel) and for a flame with L — 128 and G = 1/4 
(bottom panel). The areas below the solid and dashed lines 
are equal to hr/l- 



limit. One can understand this as follows. 

From Eq. (|6af) . we can see that vorticity is generated in 
the regions with significant temperature gradients, e.g., 
on the scale hr, and is advected by the flow on spa- 
tial scales of order hy. Thus, positive vorticity is gener- 
ated in the domains nl < x < (n + 1/2)/, while negative 
vorticity is generated in the domains (n + 1/2)1 < x < 
(n + 1)1; however, the total (signed) vorticity in the do- 
main is conserved. Diffusion of vorticity occurs predom- 
inantly across the boundaries x — nl/2. More directly, 
it is straightforward to integrate the vorticity equation 
(Ea. lfja|) over the area nl < x < (n + 1/2)1 to obtain the 
vorticity balance, 



gs 2 - 



poo 


'd 2 v y - 


+ 


'# 2 V 




' — oc 




x=0 


dx 2 





dy. 



£=1/2 



Here Cl is the total vorticity generated in the roll nl < 
x < (n + 1/2)1, and diffused through its boundaries. 
In Fig. [§] we have plotted the non-dimensional vortic- 
ity generation, averaged in the area element (1/2, Ay), 
fuj = g{l/ s l){A(l/Ay), and corresponding fluxes across 
the roll boundaries. Note that only diffusion can lead 
to vorticity transport across roll boundaries because the 
transverse flow vanishes identically at the separatrices. 

In other words, in the thin flame limit, vorticity gen- 
eration depends on the LG product, but not on the vis- 
cosity; however, in steady state, we know that vorticity 
generation and destruction must balance exactly. Since 
the vorticity destruction depends on the diffusion term 
PrV 2 w, which decreases as L increases, balance can only 
be achieved if the length of the vorticity diffusion region 



(i.e. the separatrix separating adjacent rolls) lengthens. 
Thus, we expect hy to scale inversely with v = Pr/L. 
Indeed, we expect hy — > oo as Pr — * 0. 



C. Comparison with linear stability analyses 

A thorough analysis of the linear behavior of our sys- 
tem was presented by Zeldovich et al. [2l|; in this sub- 
section, we compare our results with theirs. 

The simplest case studied is the so-called Landau- 
Darrieus instability, in which the flame is considered as 
a simple gas dynamic discontinuity [22l |23| . The fluid 
on either side of the discontinuity is assumed to obey the 
Euler equation; the fluid is assumed to be incompress- 
ible; there is no temperature evolution equation; and the 
front is assumed to move normal to itself with a given 
laminar speed. The important parameter is the degree 
of thermal expansion, 6 = /9f ue i/Pash, across the flame 
front. The resulting instability growth rate is propor- 
tional to the product of the laminar flame speed and the 
wavenumber of the front perturbation, with a coefficient 
of proportionality depending on 9. For 9 = 1, which 
corresponds to the Boussinesq limit, the growth rate is 
identically equal to zero. 

The Landau-Darrieus model is however not valid for 
wavelengths short compared to the flame thickness, for 
which it predicts the large st g rowth rate; this deficiency 
was resolved by Markstein |24| , who introduced an empir- 
ical "curvature correction" for the flame speed within the 
context of the Landau-Darrieus model. One consequence 
of this correction is that the instability is suppressed for 
wavelengths shorter than a specific critical cutoff wave- 
length, while for wavelengths much larger than this cutoff 
lengthscale the growth rate approaches zero as 1/L, just 
as in the Landau-Darrieus model. 

Gravity can be introduced in this type of model in 
a very similar way, as shown by Zeldovich et al. |'2l| . 
Rewriting their result in our notation, and taking into 
account 6—1 and Le = 1 (which leads to the Markstein 
curvature correction constant being set equal to unity), 
we can reduce their final result to the following expression 
for the growth rate, 



~1 = —,k Jl + k(k-2) + 2G/k-l-k 
25 L 



(11) 



where k = 2n/L. A more elaborate model for the flame, 
introduced by Pelce & Clavin 2^, avoids the empiri- 
cal curvature correction constant and, in the Boussinesq 
limit, gives the growth rate expression 



25 



s/T+2G/k-l 



k 



v/l + 2G/fc 



(12) 



In the limit of thin fronts, L 1, both models reduce to 
the same expression, which also recovers the LG similar- 
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FIG. 10: Growth exponent for a single mode and two initial 
amplitudes A = a/8 = 4 and A = 8. The dimensionless 
wavelengths are L = 1/6 = 512 (dashed) and L = 1024 (solid). 
The dotted lines correspond to the linear stability analysis 
prediction, Eg. 1131 . 



ity scaling already discussed above, 




(13) 



To compare our calculations with this result, we have 
computed the growth rate for a single wavelength for 
a system with L = 512 and L = 1024 (see Fig. ITU)! , 
The growth rates predicted by Eq. (| 1 31) are shown as hor- 
izontal dotted lines for each LG product. An ideal sys- 
tem in the linear regime would have a constant growth 
rate; in our simulations we observe an essentially time- 
independent growth rate only after some transitional pe- 
riod, t < 0.1 1 / ' s Ql and before the flame stabilization time, 
which depends on parameters. The transitional period at 
the beginning of our simulations can be explained by arti- 
ficial initial conditions, e.g. zero velocity and prescribed 
temperature profile across interface. The decrease in the 
growth rate at later times is related to the stabilization 
of the flame front. Naturally, the faster-growing insta- 
bilities with higher LG product and the systems starting 
with larger initial amplitudes reach the steady-state more 
quickly. In addition, we observe the influence of the finite 
flame thickness — plots with L = 1024 approache closer 
to the infinitely thin limit than plots with L = 512. But 
in spite of the finite flame thickness and non-zero viscos- 
ity, one can clearly see the similarity scaling on LG and 
good agreement with theory (Fig. Illfl . 

In order to obtain the stability condition, we set 7 = 
in the expressions Ijlll) and l|12fl . and obtain 



for the Markstein model, and 



G - k 
G CI - - 



~ (l + k+ V(l + fc) 2 +4fc) 2 - 1 



(14) 



(15) 



FIG. 11: Growth exponent for a single mode measured at the 
maximum for A = 4 and L — 1024. The solid line corresponds 
to the linear stability analysis prediction, Ea. Hl.'ifl . 



for the Pelce and Clavin model. We emphasize that both 
of these models assume an inviscid fluid, while viscosity 
is present in our simulations. In Fig. I12l critical gravities 
derived using both models are plotted next to numerical 
simulations data for different Prandtl numbers. 

Similarly, we can consider the relation between the 
instability growth rate and the amplitude of the stable 
flame front, using the assumption that the flame front is 
composed of joined parabolic segments whose amplitude 
is small when compared to their wavelength [21] ]. The 
resulting estimate depends on the growth rate, Eq. i|13fl , 



1 



I 



I 8 Vs. 



Comparing the result with the fit derived from the ex- 
perimental data shown in Fig. [3] we notice that, in the 
thin front limit and for values of G larger than critical, 
both numerical experiment and theoretical model predict 
h T /l « 0.22 Vl~G. 

Finally, we note that a quick comparison of the 
asymptotic behavior of the Rayleigh- Taylor 0, 0] and 
Landau-Darrieus [22I |23| instabilities for large L gives 
7 oc L~ x / 2 for Rayleigh- Taylor and 7 oc L^ 1 for Landau- 
Darrieus. In our Boussinesq case, the same asymp- 
totic limit gives 7 oc L~ 1//2 : the instability behaves like 
the Rayleigh- Taylor instability at long wavelengths, i.e., 
longer wavelengths grow more slowly, but saturate later 
and reach larger front speeds. 



D. Transition to the travelling wave 

The transition time during which the temperature 
front is formed is of the order of the laminar burning 
time across the period, 0.5 l/s 0l but a much longer time is 
needed to stabilize the velocity pattern behind this front. 
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FIG. 12: Critical gravity G cr for different values of Pr = 
v/k. Our results are fit with power laws, in the form G cr = 
C(2iv/L) n , with measured values of n = 2.68, C = 43.94 for 
Pr = 4; n = 2.01, C = 10.31 for Pr = 1; and n = 1.73, 
C — 6.21 for Pr = 1/4. The two solid lines are provided by 
inviscid theory (Pr = 0), corresponding to the Markstein and 
the Pelce & Clavin models. 



Fig. 03 illustrates the process for a moderate value of L; 
in Fig. 1131 we show snapshots for a flame with an L value 
closer to the Ray leigh- Taylor limit just discussed. In- 
deed, Fig. 1131 shows morphology strongly reminiscent of 
the Rayleigh- Taylor instability, namely upward-moving 
"bubbles" and downward moving "spikes" . As mentioned 
earlier, the reaction stabilizes the shape of the moving 
front, and eventually the flame interface will become 
smooth, similar to those shown in Fig. [21 The typical 
flame stabilization time is of the order of L/s a , provided 
the initial perturbation amplitude is large enough (on 
the order of a fraction of L). During the transition, the 
system with larger L/5 ratio develops more complicated 
structures (compare Fig.^Jwith Fig. [8} — the details on 
the scale of flame thickness and smaller are consumed by 
the burning. A similar effect is observed in the Rayleigh- 
Taylor instability on the dissipation scale, but due to 
viscosity and diffusion rather than burning. 

The images shown in Fig. 1141 illustrate the propagation 
of a flame with eight wavelengths (with L = 16, G = 4) 
within the computational box with reflecting boundary 
conditions. The chosen parameters place the system well 
inside the unstable regime, and, by the time t ps 30<5/s o 
the system forms the curved travelling wave solution with 
wavelength L — 16. This solution is exactly the same as 
the curved solution obtained in the half-wavelength com- 
putational box, propagates with the same speed, and re- 
mains unchanged until time t « 100 S/s - (We note that 
the "wall effect" seen in this figure reflects both the pres- 
ence of the walls (and choice of boundary condition at 
the walls) and the choice of phase for the initial pertur- 
bation.) 

The symmetry of the initial conditions requires zero 



horizontal velocity at x = nl/2, n = 0, 1, 2, this sym- 
metry constraint is clearly broken for t > 1005 /s , and 
the travelling wave solution becomes violently unstable. 
The cause of this symmetry breaking is apparently accu- 
mulated numerical errors (noise) in the calculation. The 
source of this noise is the iterative solution for the stream 
function, so the noise has the wavelength of the compu- 
tational domain. Since perturbations with larger wave- 
lengths move faster, the system will eventually pick the 
travelling wave configuration corresponding to the largest 
possible wavelength — in the example shown, the wave- 
length L = 256 (twice the box size). 

The instability shown in Fig. ^] is not related to 
metastable behavior near G cr discussed in [T^j — both 
wavelengths present in the system are unstable for G = 4. 
Rather, this simulation is an illustration of the fact that 
the small wavelengths have faster initial growth rates, 
but saturate at lower speeds. As a result, the instability 
exhibits a strong inverse cascade. More careful modelling 
of the noise introduced to the system, as well as more re- 
alistic treatment of boundary conditions at the walls, will 
be necessary to learn about the instability dynamics in 
an infinitely large domain; in particular, we believe that 
periodic boundary conditions should be imposed at the 
walls in order to study this problem further. In such a 
system we would expect unbounded growth of instabil- 
ity size; in a natural system we would expect the upper 




FIG. 13: The temperature, vorticity generation rate, and 
vorticity (from left to right) for the system with L — 256 
and G = 4 at time t = 72 5 /s a . The initial amplitude of the 
perturbation was a/ 1 = 1/8. 
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FIG. 14: Symmetry breaking due to numerical noise and resulting instability. The snapshots are taken at times ts /5 = 
100, 110, 120, 130, 140, 150. 



bound to be set by extrinsic spatial scales of the physical 
system. 



IV. SUMMARY AND DISCUSSION 



stability we observe seems to involve seeding, and strong 
growth, of modes with wavelengths much larger than the 
wavelength of the dominant front disturbance. We are 
currently investigating this instability in greater detail. 



In this paper, we have studied the fully nonlinear be- 
havior of diffusive pre-mixed flames in a gravitationally 
stratified medium, subject to the Boussinesq approxima- 
tion. Our aim was both to compare our results for a 
viscous system with analytical (and empirical) results in 
the extant literature, and to better understand the phe- 
nomenology of fully nonlinear flames subject to gravity. 

The essence of our results is that the numerics by and 
large confirm the Markstein and Pelce & Clavin mod- 
els, and extend their results to finite viscosity. We have 
shown explicitly that there is an extended regime for 
flames with finite flame front thickness for which the scal- 
ing on the LG product applies (as it is known to do in 
the thin flame front limit). We have also examined the 
details of the flame front structure, and are able to give 
physically-motivated explanations for the observed scal- 
ings, for example, of the flow length scale behind the 
flame front on Prandtl number. 

We have also observed a potentially new instability, 
which arises when noise breaks the symmetry constraint 
of the initial front perturbation. Our study suggest that 
this instability differs significantly from finger merging 
behavior of the non-linear Ray leigh- Taylor instability, in 
which the finger merging process resembles a continuous 
period-doubling phenomenon (e.g. adjacent fingers at 
any given generation merge in pairs). In contrast, the in- 



Finally, it is of some interest to consider the implica- 
tion of our results for astrophysical nuclear flames, as 
arise in the context of white dwarf explosion. Using the 
results of Timmes & Woosly |2(j , we find that we would 
be far into the thin flame limit, with a density jump 
at the flame front Ap ~ 0.1/9; hence our Boussinesq re- 
sults are rather marginal in their applicability. Never- 
theless, one can ask what the expected flame speed up 
would be in this limit; using our results we find that 
s/s re (1 + 0.0486 LG) 1 / 2 , with LG = (Ap/p)lg/s 2 . Us- 
ing the length scale of the order of a fraction of white 
the dwarf radius, I pa 10 3 fcm, gravitational accelera- 
tion on the surface of the star, q pa 10 3 fcm/s 2 , and 
the laminar flame speed given by [2q . s ~ 100 km/s, 
we obtain LG ps 10, and consequently a speed up of 
s/s pa 1.2. Smaller laminar flame speeds would lead 
to the flame velocities independent of the laminar flame 
speed, s = 0.23 (IgAp/p) 1 / 2 ps 100 km/s, which could 
be also derived using the rising bubble model [3. Ev- 
idently, the flame speedup in this limit is very modest. 
Whether compressibility has much effect on this conclu- 
sion remains to be established and is now under active 
investigation. 
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